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Abstract 

We propose the use of mixing strategies to accelerate the convergence of the common iterative algorithms 
utilized in Quantum Optimal Control Theory (QOCT). We show how the non-linear equations of QOCT can 
be viewed as a "fixed-point" non-hnear problem. The iterative algorithms for this class of problems may 
benefit from mixing strategies, as it happens, e.g. in the quest for the ground-state density in Kohn-Sham 
density functional theory. We demonstrate, with some numerical examples, how the same mixing schemes 
utilized in this latter non-linear problem, may significantly accelerate the QOCT iterative procedures. 
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I. INTRODUCTION 



Quantum optimal control theoryUJ, |2i |3i |4|] (QOCT) answers the following question: A sys- 
tem can be driven, during some time interval, by one or various external fields whose temporal 
dependence is determined by a set of "control" functions. Given an objective (e.g., to maximize 
the transition probability to a prescribed final state, the so-called target state), what are the control 
functions that best achieve this objective? 

In the more general context of dynamical systems, optimal control theory is widely used for 
engineering problems, and its modem formulation was established in the 1950's.|| 5|] Th e trans- 
lation of these ideas to Quantum Mechanics was initiated in the 1980's.|^ R Isl. IX 10. [ill [l2[ 
[isll Recently, the field has received increasing attention due to the parallel 
advances in experimental control techniques: femto - and atto - second laser sources with pulse 
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shaping. 11191 120L 12111 and learning loop algorithms. II22I1 These new developments call for corre- 
sponding theoretical efforts. 

The computational solution of the QOCT equations may impose an enormous burden. Any 
algorithm requires multiple forward and backward propagations of the quantum system under 
study. This can be very cumbersome, depending on the level of theory employed to model the 



process. The development of efficient a 



ficient schemes already exist. [1231 
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gorithms is therefore essential. And, in fact, rather ef- 



26D The most effective choices are closely related and 



can be grouped in a unified framework. [|27ll The equations to be solved are non-linearly coupled 
initial-value partial differential equations, and must be solved iteratively. These iterative proce- 
dures can be described in the following way: one input field is passed to an "iteration functional" 
that tests its performances and produces an improved "output" field. This output field can then be 
used as input for the iteration functional. Upon solution, output and input fields coincide at the 
"fixed-point" of the iteration functional. 

We must therefore search for the fixed-point of some non-linear functional. One prominent 
example of this kind of fixed-point problems is the Kohn-Sham (KS) formulation of density- 
functional theory (DFT).[28] In this field, it was soon realized that the naive use of the output 
produced in one iteration as input for the next one leads to poor (or no) convergence, and this 



observation suggested the use and development of "mixing" techniques: [|29L 30 . 3ll] the input for 
each iteration is a smart combination of the output of the previous iteration and several inputs or 
outputs of former iterations. The result is typically a very significant acceleration in the conver- 
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gence - and even the possibility of finding a solution in cases where no mixing (or trivial "linear" 
mixing) is unable of finding one. 

In this work, we propose the use of those mixing strategies to accelerate the convergence of the 
iterative algorithms used in QOCT. We demonstrate how they can significantly reduce the iteration 
count - yet the performance and degree of gain, of course, depends on the details of each particular 
model. The procedure should be viewed as a scheme to accelerate (and not substitute) the existent 
iterative algorithms; in particular, it will be made evident that the mixing should be switched on 
after a couple iterations have been made and the control function is not too far away from the 
solution - fortunately, it is precisely the regime where the existent algorithms behave better. 

The description of the proposed methodology is provided in Section [III Some numerical evi- 
dence supporting the advantages of its use is shown in Section [nil Atomic units are used through- 
out. 

11, METHODOLOGY 

We recall the essential equations of QOCT, making no attempt to state them in full generality - 
the basic ideas can be generalized in different ways suitable for a broad class of situations; however 
the reader should find no difficulties to translate our suggested enhancements to those variations. 

We consider a system characterized in the absence of external fields by a Hamiltonian ^o- One 
external "control" operator t{t)V may drive it during some time interval, where e(f) is a "control 
function". The system is therefore governed by: 

H{t)=HQ + z{t)V , 0<t<T, (1) 

which drives the system from its initial state |^o) to a final state |^(r)). The purpose of the 
QOCT algorithms is to find that e{t) that maximizes the value of a target - in mathematical terms, 
a functional of the evolution of the state, 7i [^] . hi many cases it depends only on the value of the 
state at the end of the propagation. And in most cases it takes the form of the expectation value of 
some operator O: 

7i[^] = (^(r)|d|^(r)). (2) 

In order to produce a physically meaningful process, the maximization of 7i must be con- 
strained: on one hand, one should limit the search space of £; on the other hand one must ensure 
that the evolution \^{t)) indeed follows from Schrodinger's equation. In mathematical terms, this 
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translates into the maximization of the functional: 

7[^,X,£] = +72[e] +73[^,X,e] , (3) 

rT 

72[e] = -a / d?£2(0, (4) 

JO 

73[^,X,£] = -2JmyJ^i(x(0|i^-^o-e(OV|^(0)- (5) 

The 72 functional penalizes the "fluence" of the control field - ensuring that the maximization pro- 
cedure does not lead to infinite values. The 73 functional ensures that ^{t) satisfies the Schrodinger 
equation, and it introduces a new "Lagrange-multiplier" wave function, %. The Euler-Lagrange 
equations satisfied at the stationary points of 7 are: 

(6) 
(7) 
(8) 
(9) 
(10) 

Numerous modifications and extensions to these equations are possible; for example, the possibil- 
ity to include dissipation. B2il to account for multiple objectives. [|33ll to deal with time-dependent 
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= Jm(%(0|V|^(0)- 



targets, yj. 
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to add spectral and fluence constrains. 113711 or to work with more general inho- 
mogeneous integrodifferential equations of motion. ||38|1 In order to keep the present discussion as 
simple as possible, we have chosen to present this "standard" set of equations, but we stress that 
the algorithmic enhancements discussed below may be applied to the modified versions. 

The control equations (l6l)-(fT0l) are coupled and one must look for a self-consistent solution. 
This requires an iterative scheme, which we choose to write here, for reasons that will become 
clear later, in "functional" form: an "iterator" F takes an input control function £ and produces an 
output, which is then used as input for the following iteration. 

The simplest option would be what can be called a "straight iteration": given a trial control 
function £'^*) {k is the iteration index), the output F[£(^)] is constructed by taking the steps: 

1. Propagate, from |^(0)) = |^o) to |^(r)) with £W. 

2. Propagate backwards, from \x{T)) = d\^'{T)) to |%(0)), also with £^^\t). During the evo- 
lution, calculate the output field F[e^^^]: 

aF[£«](0=Jm(%(0|y|^W)- (11) 
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3. Define e^^+i) = F[e^^^, and repeat from step 1 until convergence is reached (F[e] = £). 



This procedure was already used, for example, in the seminal work of Kosloff et al UM- As 



discussed by Somloi et al [|26|], doing such a straight iteration in general does not lead to conver- 
gence. One possible way to cure this problem is to set e^^+i) = £(^) + yF[£^^^]; the parameter y may 
be set by performing a line-search optimization, such that e'-'^+^'' produces the maximal objective 
J. This idea is in fact a first-order approach to the schemes discussed below. 



One monotonically convergent algorithm was introduced in Ref. [|24ll - we will refer to this 
algorithm as ZR98. It can be described in the following way: given the trial control function e^^\ 
the output F[£(^)] is constructed by taking the steps: 

1. Propagate, from |^(0)) = |^o) to with E^'^). 

2. Propagate backwards, from |%(r)) = d\'^{T)) to |x(0)), with £ defined as: 

ae(0 = 3fm(x(0|t>|^(0). (12) 

E must be obtained "on the fly", from the values of the propagating \x{t)) and the previously 
obtained \^{t)). 

3. Propagate forward, from |^'(0)) = |^o) to |^'(r)), using the output field F[£(^)] (t), which 
is now defined as: 

aF[£(^)](0 = 2rm(x(0|V|^'W)- (13) 

One can then simply define £(*+i) = F[£(*)], and proceed to the next iteration (note that in this 
case one does not need to perform explicitly step 1 again, since it repeats step 3 in the previous 
iteration). The solution, i.e. the optimal field, is obtained when the iteration finds a fixed point: 

F[£] =£. 

At this point, some remarks are in order: 



The seminal algorithmic work presented in Ref. [ |23l ] deals with a slightly modified version 
of the previous scheme, suitable for a particular (but very common) type of target - the 
objective is to maximize the population of a given target state ^target, and therefore the 
operator O is the projection onto that state. In this case, the generating functionals can be 
defined in such a way that the propagations for ^ and % are decoupled (note that in Eqs.lMTOl 
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^ and % are coupled through Eq. |9l): 



at 


[Ho + £{t)V] 1^(0), 


MO)) = 


l^o), 




[Ho + Eit)V] IxW), 


\X{T)) = 


^target) , 


a£{t) = 


Jm(x(0|V|^(0)- 



(14) 
(15) 
(16) 
(17) 
(18) 



This algorithm, presented in Ref.|23|, is essentially the same as the one defined by Eqs.lMTOl 
except for the fact that the initial value point for the propagation of % is now given by Eq.fTTi 
In the following, we will refer to this scheme as ZBR98. 



can also be viewed as a 



271 both schemes can be 



• The formulation of Krotov as implemented by Tannor et aln25 
modification of the previous scheme and, as discussed in Ref 
written as members of the same family. 

• A remarkable property of both ZR98, ZBR98 and Krotov's schemes is their monotonic con- 
vergence. Also, and most importantly, they typically provide very fast improvements during 
the first iterations - when the trial field is very far from the solution field. Unfortunately, 
this rate of convergence slows down as the iteration count grows. 

• The previous description of the Z(B)R98 algorithms suggests that the cost, per iteration, is 
that of two wave function propagations. However, this can only be possible if, at every time 
step in either the forward or backward step, the wave functions are stored, and then retrieved 
from memory when performing, respectively, the following backward or forward step. In 
practice, this can be time consuming, and the best way is actually to propagate once again 
with the same field, but in the opposite time direction. In this case, one needs, in fact, four 
propagations per iteration. 

We propose now to utilize the simpler straight iteration, but with an important change: We do 
not use e(^+i) =F[eW]. Instead, we can use the sophisticated mixing algorithms that have proved 
so useful in the field of Kohn-Sham density functional theory. [|28|l such as, for example, the one 
presented in Ref. OOll (which we will call modified Broyden's algorithm). Other options would 
be equally valid - for example the work presented in Ref. [13 ill . One can use them exactly in the 
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same way as they are used in electronic structure calculations. The only external ingredient that 
the algorithms necessitate, and which is different for each problem, is a dot product definition for 
the relevant variables, which in the QOCT case are the control functions. We take the obvious 
choice: 

(£l|£2)= / dt£i{t)e2{t) . (19) 

Jo 

In essence, the gist of Broyden's scheme (and of all of the other so-called "mixing" strategies) 
consists of making use, in order to define E^^^^-', not only of F[£^^^], but also of a number s of 
previous iteration values: 

e^'+^) = Gn^xing[ {e^'-^\F[£^'-j^]yj-J,]. (20) 

The functional Gmixing is chosen in some way designed to minimize the distance D between input 
and output: 

D{F[e],£) = {F[e]-e \ F[£] -e)^!/^) . (2i) 

These functionals are essentially based on approximations to the conventional Newton- 
Raphson iteration. Let us define T[e] = F[e] — £; in the vicinity of E^^) (the k-th iteration ap- 
proximation to the solution), the functional T can be linearized: 

r[E]^r[E(^)]+7W[E-E(^)], (22) 

where is the Jacobian of T evaluated at E*^^-*. This can be rewritten as: 

e_e«+G(^)[r[E]-r[E(^)]] =0, (23) 

where G^^^ = —{J^^^)^^. Newton's iteration follows inmediately from this formula by assuming 
e(^+i) to be the solution vector (r[£(*+i)] = 0): 

Since the Jacobian (let alone its inverse) may be difficult to compute, quasi Newton-Raphson 
schemes utilize approximations to it; in Broyden's family of schemes, these are also built itera- 
tively. Therefore, the matrices G*^^-* do not verify Eq. [23] until convergence. 



Johnson's proposal, yOD in particular, consists of generating G*^^+^) by minimizing the follow- 
ing functional: 

k 

E = a)^||G(*+i) -G^ll + £ co^IAeW + G(^+i)ArW|2, (25) 



where Ae^ = £("+!)-£(«), AjW = T[e^"+'^^] - r[£W], and (O,; are a set of real positive constant 
weights. The idea is therefore to minimize the error in the inverse Jacobian (first term in the 
definition of E), at the same time making sure that the new guess for G verifies Eq.|23]as closely 
as possible - not only for the last iteration E^^), but also for all the previous ones. 

The minimization of E with respect to G^^+^^ leads to an iterative formula; this formula can 
then be plugged into Eq. [241 The final result reads: 

e(^+i) = + G(i)r[e«] - £ (DnjknU^"^ , (26) 

where: 

= G(i)ArW+A£W, (27) 

Y«=£«.(ArW|rW)|3„,, (28) 

n=l 

%i = {(oll + a);J, (29) 

aij = (Oi(Oj{AT^j'^\AT^'^). (30) 

In this manner, the new field £(^+i) can be obtained from information gathered in the previous 
iterations; it only remains to choose an appropriate initial guess for the inverse Jacobian, G*^^h we 
merely set it equal to some constant times the identity, a/. This constant a can be freely chosen 
(it ultimately determines the amount of "output" field to be utilized in the mixture, as it can be 
understood inspecting Eq. [26l) . and it is a matter of experience to determine a reasonable value - 
likewise for the (O,- constants. In typical DFT codes, the CO, constants are not adjusted for each run; 
the same values are used for all systems, a is usually set initially to some "aggressive" value (i.e. 
large value, meaning a large proportion of the "output" is used), and, if convergence is not found, 
it is reduced in a subsequent run. We expect that the same strategy should hold for QOCT runs. 

A careful analysis of Eq. [26lalso reveals that we do not need to manipulate or store objects of 
size A^^, where A'^ is the dimension of the problem field £. The cost of the operations is of order ^A^^, 
where s is the number of previous iterations to be considered in the formula. This is the key reason 
to utilize this modification of Broyden's scheme, since in a typical QOCT problem the dimension 
A'^ is given by the number of time steps in the propagation, which can be easily of the order of 10^. 

Regarding numerical details, we have implemented the QOCT machinery in our electronic- 
structure octopus code. BQII Since one of the tasks of this code is to solve the KS DFT problem, the 
mixing strategies cited above are implemented. This platform is specialized in the time-dependent 
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version of DFT, TDDFT, and therefore contains sophisticated time-propagation schemes . 1401] uti- 
lized for the results shown below. We have recently employed our QOCT machinery - without 
making use of the mixing strategies - to model the control of electrons trapped in two dimensional 



semiconductor quantum nanostructures. ll4lll 



III. RESULTS 

A. Asymmetric double well 

As a first example, we will use a simple, but prototypical example: the transfer of a wave packet 
from one to another well in an asymmetric double well potential. The field- free Hamiltonian is: 

1 3^ -t" ^ 

The potential is depicted on top of Fig. [TJ Qualitatively similar potentials appear in many ar- 
eas in Physics, and in Quantum Chemistry they are sometimes used to model isomerization 



problems. [|42l] The external control couples to the system through the dipole operator: V = —x. 

The ground state, is localized in the left, deeper, well. The first excited state, which will be 
our target, ^target^ is localized in the right well. These states are also depicted in Fig.[Tl The target 
operator O is in this case the projection [%arget) (^target I > and therefore the preferred algorithm 
to start with is the one presented in Ref. [|23|] . ZBR98. In addition, we will use straight iteration 
assisted by Broyden's mixing. 

Fig. [2] displays four QOCT runs, each of them considering a different initial guess for the 
solution field: 

e^^\t) = Eocos{iot) . (32) 

We try the optimization with four different values of Eq, as displayed on each of the panels. Both 
the values of the total functional J and of the objective 7i are displayed. Before describing the 
results, we should point out that, as discussed above, there are some adjustable parameters to 
completely define the mixing algorithm: (i) the number s of previous iterations considered in the 
construction of the algorithm - which is set to four in this case; (ii) a number a that specifies the 
amount of output field, aF[£(^)], that is utilized in the mixing - we use a = 0.1; (iii) also, it is 
sometimes advisable to stop the algorithm every given number of steps, and restart erasing the 
memory from previous iterations - in Fig. [2l however, we have chosen to put the straightforward 
algorithm. We have made no attempt to optimize the method by taking advantange of this freedom. 



The results are very promising: except for the case (top left) where the initial laser field has 
very low amplitude (Eq = 0.01), leading to a very small initial overlap Ji, in all the other cases 
Broyden's mixing converges faster - and in fact converges to a different, better maximum. Unfor- 
tunately, the exception in the top left comer is very disappointing since the procedure yields the 
zero field - which is also a solution of the QOCT equations, but certainly not the desired one. 

And that exception is indeed specially important, since it exemplifies an important weakness 
of using straight iteration together with Broyden's algorithm: the algorithm behaves very poorly 
if the initial guess is not good enough. Fortunately, this is precisely the regime where most of the 
already existent algorithms behave better - and therefore one can devise a "hybrid" procedure: a 
few iterations with, for example, ZBR98, followed by the mixing iterations. Fig. [3] displays the 
results obtained in this way: at iteration number three, the ZBR98 procedure is stopped. Then, 
after a couple of irregular iterations (that irregularity can however be controlled by a more careful 
selection of the parameter a mentioned above), the results are significantly better. 



B. Morse potential 

For our second example, we have chosen a case that has already been discussed in the literature: 
the vibrational excitations in a Morse potential model for the OH bond. This potential function is 



given by: [|43ll 



Vix) = Do [exp(-p(.^ - ro)) - 1]' - Dq , (33) 

where, for the OH case, the parameters are chosen to be: Dq = 0.\ 994, (3 = 1.189, ro= 1.821. The 
coupling to the external function is given now by a dipole potential operator approximated by the 
function: 

V{x)=^oxexp{-x/r*), (34) 

where fiQ = 3.088 and r* = 0.6. The objective is now to populate the first excited state, starting 
from the Morse ground state. The total propagation time is T = 30000 a.u.(^ 0.725 ps); the 
initial trial input field is the zero field, and the penalty factor is a = 1 . This is precisely the first 



example discussed by Zhu et al [|23|] to demonstrate the performance of the ZBR98 algorithm. We 
have replicated those calculations with our codes, and in the following we demonstrate how the 
addition of mixing strategies significantly boosts the performance of the original scheme. [44] 

The results are shown in Fig. |4l First of all, we should note that the attempt to apply a straight 
iteration scheme right from the zero-th iteration - whether or not assisted by mixing techniques - 
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will fail, since the initial input field is already a solution to the QOCT equations. This solution, 
however, is an unstable point, and the ZBR98 in this case relies on numerical error to abandon 
this unstable solution, and then proceeds until convergence into a maximum. The behavior of this 
algorithm is summarized by the circles in Fig. HI which show both the values of Ji and J at each 
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Also, the 



iteration step. These results are almost exactly the same as the ones given in Ref. 
final converged field (shown in the inset of Fig. HJ coincides with the one reported in that work. 

In order to speed up the convergence, we utilized the modified Broyden's mixing algorithm 
to accelerate the straight iteration scheme, starting from the field obtained after the first ZBR98 
iteration. The results are shown in Fig. |4] with squares; the thick black curve in particular refers to 
the convergence of the J functional. It may be seen how the final converged value, J{°°) = 0.885, is 
obtained in a few iterations. In order to achieve the same level of convergence, ZBR98 necessitates 
around 50 iterations. Once again, it should be noted that in a usual implementation, each straight 
iteration step will be half as costly as a ZBR98 step. 

IV. CONCLUSIONS 

To summarize, ideas borrowed from a particular field of computational physics (e.g., density 
functional theory techniques) have been used successfully in the completely different context of 
QOCT algorithms. The reason for success is, of course, the underlying parallelism in the equa- 
tions, if regarded with the appropriate "abstract eye". 

This work by no means proves, in mathematical rigour, the universal advantage of using mixing 
strategies for all QOCT problems. However, we have observed important speedups in most cases 
(as in the ones presented in this article), and we feel that the analogy with the KS DFT problem 
provides convincing evidence about the usefulness of employing mixing for QOCT. Two important 
features of previous algorithms are, unfortunately, lost: the monotonic convergence (thi s could be 



cured, nevertheless, by adapting the mixing scheme proposed by Bowler and Gillan plU ). and the 
usual large gains during the first iterations when the initial guess is far from the solution. 

This excellent performance of both Z(B)R98 or Krotov's algorithm for the first iterations sug- 
gests the use of the mixing schemes not as a substitute, but as a complement - as demonstrated in 
our sample runs. Moreover, we should note that the idea of mixing several iterative steps accord- 
ing to, e.g., Broyden's scheme can also be applied on top o/the usual ZBR98, ZR98 or Krotov's 
algorithms (and not on top of the straight iteration scheme, as it has been presented here). Our 
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experience shows that, in most cases, doing this accelerates the convergence - although at the cost 
of loosing the predictable, regular, monotonic behaviour. Finally, we remark that the suggested 
algorithms can be easily mounted on top of the already working programs. 
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FIG. 1 : Asymmetric potential well (thick line), together with the initial (dotted line) and target (dashed 
Une) states. 
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FIG. 2: Convergence histories for both the ZBR98 algorithm and the straight iteration scheme assisted with 
the modified Broyden mixing scheme. Each panel displays the results obtained with a different initial guess 
(see text). 
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FIG. 3: Convergence histories for both the ZBR98 algorithm and the straight iteration scheme assisted with 
the modified Broyden mixing scheme. The modified Broyden scheme, however, is only applied after the 
third iteration. 
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FIG. 4: Convergence histories for both the ZBR98 algorithm (Unes with circles) and for the straight iteration 
scheme assisted with the modified Broyden mixing scheme (lines with squares), for the case of the Morse 
potential. Both the values for the Ji ("target functional", lower lines) and for the total J functional (upper 
lines) are shown. Inset: optimized control field. 
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